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WAVELET  SPECTRAL  EINITE  ELEMENTS  FOR  WAVE 
PROPAGATION  IN  COMPOSITE  PLATES  WITH  DAMAGES 


S.  Gopalakrishnan,  Indian  Institute  of  Science,  Bangalore,  India 

ABSTARCT 

For  the  Current  Year  (2009-2010),  A  new  spectral  plate  element  (SPE)  for  laminated 
composite  using  Daubechies  compactly  supported  wavelet  basis  functions  is  developed 
to  analyze  wave  propagation  in  an  anisotropic  laminated  composite  media.  The  element 
is  based  on  the  Classical  Laminated  Plate  Theory  (CLPT).  The  element  is  formulated 
using  the  recently  developed  methodology  of  spectral  finite  element  formulation  based  on 
the  solution  of  a  Polynomial  Eigenvalue  Problem  (PEP).  By  virtue  of  its  frequency- 
wavenumber  domain  formulation,  single  element  is  sufficient  to  model  large  structures, 
where  conventional  finite  element  method  will  incur  heavy  cost  of  computation.  First  the 
Wavelet  spectral  element  procedure  is  outlined,  which  is  followed  by  a  section  on 
wavenumber  computation  as  a  function  of  fiber  directions  is  given.  A  number  of 
numerical  examples  are  provided  to  show  the  efficiency  of  the  formulated  element.  In 
most  cases,  the  results  from  the  wavelet  formulations  is  compared  with  the  conventional 
finite  elements  to  show  the  computational  superiority  of  the  formulated  wavelet  spectral 
element  for  certain  class  of  problems. 

1.  INTRODUCTION 

Composites  are  being  used  in  aircraft  structures  because  of  their  several  favorable 
properties.  Study  of  wave  propagation  in  such  composite  structures  is  of  much  relevance 
because  it  helps  one  to  understand  their  behavior  under  high-frequency  impact  load 
encountered  in  gust,  bird  hit,  tool  drop,  etc.  Apart  from  this,  wave  propagation  analysis 
finds  important  applications  in  structural  health  monitoring  using  diagnostic  waves, 
control  of  noise,  and  vibration.  However,  the  behavior  of  composites  at  high  frequencies 
is  more  complicated  than  it  is  for  their  metallic  counterpart  because  of  the  presence  of 
anisotropy,  and  very  little  literature  is  specifically  available  on  transient  wave 
propagation  analysis  of  composite  structures. 

Wave  propagation  deals  with  loading  of  high-frequency  content  and  FE  formulation  for 
such  problems  is  computationally  prohibitive  because  it  requires  large  system  size  to 
capture  all  the  higher  modes.  These  problems  are  usually  solved  in  the  transformed 
frequency  domain  using  Fourier  Transform  based  methods  and  FFT  Based  Spectral 
Element  formulation  (ESEE)  [1,  2]  is  one  such  method,  especially  tailored  for  wave 
propagation  analysis.  In  FSFE  formulation  for  2D  structures  [1],  nodal  displacements 
are  related  to  nodal  tractions  through  a  frequency-wave-number  dependent  stiffness 
matrix.  The  mass  distribution  is  captured  exactly,  and  the  accurate  elemental  dynamic 
stiffness  matrix  is  derived.  Consequently,  in  the  absence  of  any  discontinuities,  one 
element  is  sufficient  to  model  a  plate  structure  of  any  length,  but  unbounded  along  the 
other  lateral  direction. 


The  main  drawback  of  FSFE  is  that  it  cannot  handle  waveguides  of  short  lengths.  This  is 
because  the  required  assumption  of  periodicity  in  time  approximation  results  in  a 
''wraparound”  problem  for  a  smaller  time  window,  which  totally  distorts  the  response. 
However,  in  the  proposed  wavelet  based  Spectral  element  formulation  (WSEM) ,  the  use 
of  Daubechies  compactly  supported  wavelets  [3]  with  localized  basis  for  temporal 
approximation  removes  the  signal  wraparound  problem  and  can  efficiently  model 
undamped  finite  length  waveguides.  In  addition,  for  2D  problems,  FSFEs  [1]  are 
essentially  semi-infinite,  i.e.,  they  are  bounded  only  in  one  direction.  Thus,  the  effect  of 
one  lateral  boundary  cannot  be  captured,  and  this  can  be  attributed  to  the  global  basis 
functions  of  the  Fourier  series  approximation  of  the  spatial  dimension.  The  formulated 
2D  WSEM  also  overcomes  the  above  problem  and  can  accurately  model  2D  plate 
structures  of  finite  dimensions.  This  is  again  due  to  the  use  of  localized  Daubechies 
scaling  functions  as  the  basis  for  approximation  of  the  spatial  dimension. 

The  steps  followed  in  2D  WSFE  formulation  are  as  follows.  Here,  first  Daubechies 
scaling  functions  are  used  for  approximation  in  time  and  this  reduces  the  governing 
partial  differential  equation  (PDE)  into  a  set  of  coupled  PDEs  in  spatial  dimensions.  The 
wavelet  extrapolation  technique  [4]  is  used  for  adapting  wavelets  in  finite  domain  and 
imposition  of  initial  conditions.  The  coupled  transformed  PDEs  are  decoupled  through 
eigen  analysis.  This  temporal  approximation,  imposition  of  initial  conditions,  and 
decoupling  of  the  reduced  PDEs  are  very  similar  to  that  done  in  WSEM  formulation  for 
ID  waveguides  [5].  Next,  each  of  these  decoupled  PDEs  are  further  reduced  to  a  set  of 
coupled  ODEs  by  using  the  same  Daubechies  scaling  functions  for  approximation  of  the 
spatial  dimension.  Unlike  the  temporal  approximation,  here,  the  scaling  function 
coefficients  lying  outside  the  finite  domain  are  not  extrapolated  but  obtained  through 
periodic  extension  for  unrestrained,  i.e.,  free  lateral  edges.  Each  set  of  ODEs  are  also 
coupled,  and  decoupling  is  again  done  using  eigen  value  analysis. 

Presence  of  elastic  coupling  in  anisotropic  laminated  composite  plate  results  in  coupled 
governing  differential  equations.  Thus,  the  final  reduced  ODEs  after  spatial  and  temporal 
approximations  are  in  the  form  of  a  set  of  coupled  ODEs,  for  each  temporal  and  spatial 
sampling  point.  The  solution  of  these  ODEs  to  derive  the  exact  shape  function  involves 
determination  of  wavenumbers  and  the  amplitude  ratio  matrix.  Unlike  isotropic  cases, 
here  the  process  of  solution  is  more  complicated  and  is  done  by  posing  it  as  polynomial 
eigenvalue  problem  (PEP).  It  should  be  mentioned  here  that  similar  to  2D  FSEM,  the 
frequency-dependent  wave  characteristics  corresponding  to  each  lateral  (T)  wavenumber, 
can  be  extracted  directly  from  the  present  2D  WSEM  formulation.  However,  unlike 
FSEM,  the  wavenumbers  will  be  accurate  only  up  to  a  certain  fraction  of  Nyquist 
frequency  [6]. 

2.  Daubechies  Compactly  Supported  Wavelets 

In  this  section,  a  concise  review  of  orthogonal  basis  of  Daubechies  wavelets  [3]  is 
provided.  Wavelets  form  a  compactly  supported  orthonormal  basis  for  L^(R).  The 

wavelets  and  associated  scaling  functions  (p^  ^.  {t)  are  obtained  by  translation  and  dilation 
of  single  functions  y/{t)  and^(0 ,  respectively. 
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The  scaling  functions  (p{t) ,  are  derived  from  the  dilation  or  scaling  equation, 
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and  the  wavelet  function  is  obtained  as 
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where,  are  the  filter  coefficients,  which  are  fixed  for  a  specific  wavelet  or  scaling 
function  basis.  For  compactly  supported  wavelets,  only  a  finite  number  of  are 
nonzero.  The  filter  coefficients  are  derived  by  imposing  certain  constraints  on  the 
scaling  functions. 
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Let  Pj(f)it)  be  the  approximation  of  a  function /(t)  in  L  (R)  using  as  the  basis, 

at  a  certain  level  (resolution)  j,  then 


^)(/)(^)  =  2  ^  e  Z  (5) 
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where  cj^k  are  the  approximation  coefficients. 

3.  Reduction  of  Wave  Equations  to  ODEs 

3.1  Governing  Differential  Equations 

Using  classical  laminated  plate  theory  (CLPT)  [8]  the  three  governing  equations  with 
respect  to  the  three  degrees  of  freedom  and  w  are  given  below  (Eqns(6-8)).  Here, 

uo(x,  y  ,  t),  Vo(x,  y  ,  t),  and  w{x,  y  ,  t)  are  the  axial  and  transverse  displacements  in  x,  y  and 
z  directions,  respectively,  along  the  mid  plane,  which  is  at  z  =  0. 


(7) 


B]]6^Uq  {B]2+ 


dx' 


dxP^  (?  V  dxd  v^ 

ly  ly 

D I  I  2(Oi2"I"  ^Oci)  ^ W  DyX^ 


+  ■ 


dx^ 


dx^  d 


'11^ 

d\^ 


/  (5^M’  <5^v  ' 

=  4'^’  -  ^21,  TT  +  Ti 
\  dx  dy 


(8) 


The  stiffness  coefficients  A;/,  By,  Dy  and  the  inertial  coefficients  Iq,  h  are  defined  as 
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where  Qy  is  the  stiffness  constant  of  the  lamina  and  p  is  the  mass  density.  The  associated 
boundary  conditions  for  edges  parallel  to  F-axis  are 


2i]|4?Hg 


dx 


dv 


dx^ 


d\ 


,2 


(9) 


_  j  if  ^  ^  'l  _ 


*'^¥-'^66 


\  (?V 


dx 


dx  d\ 


(10) 


B]idUQ  B^2^^0 , 0]]ci^irx  D]2^\[ 

M,,  =  - - - - + - ^  ^  (11) 

dx  dv  dx  dv 

b-- 

,,  ^12^t'’0  D]2^W  /2(?W 

.2  +  ,  ,  -  ,  ^  .  2  +~ - 

dx  dxdy  dX'  dxdy  dx 

where  Nx  and  Ny  are  the  normal  forces  in  x  and  y  direction,  respectively.  My  is  the 
moments  about  y-axis.  The  shear  resultant  or  the  Kirchoff  shear  [1]  Vis  obtained  as 
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dv 
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where  Q  is  the  transverse  shear  force  in  the  z  direction.  Next,  governing  PDEs  and  the 
associated  boundary  conditions  derived  here  are  reduced  to  a  set  of  ODEs  using 
Daubechies  scaling  function  approximation  in  time  and  one  spatial  (Y)  dimension. 


3.2  Temporal  Approximation. 

The  first  step  in  the  formulation  of  the  2D  WSEM  is  the  reduction  of  each  of  the  three 
governing  differential  equations  given  by  Eqs.  (6-8)  to  a  set  of  PDEs  by  Daubechies 
scaling  function-based  transformation  in  time.  The  procedure  is  exactly  similar  to  that  in 
the  formulation  of  the  ID  WSEE  [5].  However,  the  key  steps  are  stated  here  very  briefly 


for  completeness.  Let  mo(x,  j  ,  t)  be  discretized  at  n  points  in  the  chosen  time 
window[0t^  ]  Let  T  =0,1, . . . ,  n-1  be  the  sampling  points,  then 

t-  AtT  (14) 


where  At  is  the  time  interval  between  two  sampling  points.  The  function  mq  (x,  ,  t)  can 

be  approximated  by  scaling  function  (p{t)  at  an  arbitrary  scale  as 
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where  y)  (referred  as  uqu  hereafter)  are  the  approximation  coefficients  at  a  certain 

spatial  dimension  x  and  y.  The  other  displacements  vo  (x,  y  ,  t),  w{x,  y  ,  t)  can  be 
transformed  similarly.  Substituting  these  approximations  in  Eqn.  (6),  using  the 
orthogonality  property  of  the  translates  of  the  scaling  functions  and  the  definition  of 
connection  coefficients  [8],  the  transformed  coupled  PDEs  obtained  are  of  the  form 
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where  T*  is  the  first-order  connection  coefficient  matrix  obtained  after  using  the  wavelet 
extrapolation  technique  [4].  These  coupled  PDEs  are  decoupled  using  eigenvalue  analysis 
of  r*  as  given  in  Reference  [5].  The  final  decoupled  form  of  the  reduced  PDEs  given  in 
Eqn  (16)  is 
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where  u^j  and  similarly  other  transformed  displacements  are 

(IS) 

where  4>  is  the  eigenvector  matrix  of  E'  and  /j  are  the  corresponding  eigenvalues. 

Following  exactly  similar  steps,  the  two  other  governing  differential  equations  (Eqs.  (7) 
and  (8))  and  the  force  boundary  conditions  (Eqs.(9)-(12))  are  transformed  to  PDEs  in  x 
and  y. 

It  should  be  mentioned  here  that  the  sampling  rate  At  should  be  less  than  a  certain  value 
to  avoid  spurious  dispersion  in  the  simulation  using  WSEM.  In  Reference  [6],  a 
numerical  study  has  been  presented  from  which  the  required  At  can  be  determined 


depending  on  the  order  N  of  the  Daubechies  scaling  function  and  frequency  content  of  the 
load. 

3.3  Spatial  (Y)  Approximation 

As  said  in  Sec.  1,  the  next  step  involved  is  to  further  reduce  each  of  the  transformed  and 
decoupled  PDEs  given  by  Eq.  (17)  (similarly  for  the  other  transformed  governing 
differential  equations  corresponding  to  (7)  and  (8))  for  7=0,1, .  . . ,  n-1  to  a  set  of  coupled 
ODEs  using  Daubechies  scaling  function  approximation  in  one  of  the  spatial  (Y) 
direction.  Similar  to  time  approximation,  the  transformed  variable  u^j  can  be  discretized 

at  m  points  in  the  spatial  window  [0,  Ly],  where  Ly  is  the  length  in  Y  direction.  Let 
^  =  0,1, . ,m  - 1  be  the  sampling  points,  then 


y  =  A  Y( 


(19) 


where  AF  is  the  spatial  interval  between  two  sampling  points.  The  function  u^j  {x,  y)  can 
be  approximated  by  scaling  function  at  an  arbitrary  scale  as 
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where  (x,  y)  (referred  as  u^ij  hereafter)  are  the  approximation  coefficients  at  a  certain 
spatial  dimension  x.  The  other  displacements  (x, y),  (x, y)can  be  similarly 

transformed.  Eollowing  similar  steps  as  the  time  approximation,  substituting  the  above 
approximations  in  Eq.  (17)  and  taking  inner  product  on  both  sides  with  the  translates  of 

scaling  functions  where  i  =0,1,  .  ,m-l  and  using  their  orthogonal 

properties,  we  get  m  simultaneous  ODEs  as  follows 
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where  N  is  the  order  of  Daubechies  wavelet  and  ^  and  are  the  connection 

coefficients  for  first-  and  second-order  derivative  defined  in  Reference  [8] 


It  can  be  seen  from  the  ODEs  given  by  Eq.  (21),  that,  similar  to  time  approximation,  here 
also  certain  coefficients  M^,^  near  the  vicinity  of  the  boundaries  (/=0  and  i=m-l)  lie 


outside  the  spatial  window  [L]  defined  by  /=0,1,  .  .  .  ,m-l.  These  coefficients  must  be 
treated  properly  for  finite  domain  analysis.  However,  here,  unlike  time  approximation, 
these  coefficients  are  obtained  through  periodic  extension,  but  only  for  free  lateral  edges, 
while  other  boundary  conditions  may  be  imposed  quite  differently  using  a  restraint  matrix 
[9].  The  unrestrained,  i.e.,  free-free  boundary  conditions  may  also  be  imposed  in  a  similar 
way  using  a  restraint  matrix,  but  interestingly,  it  has  been  seen  from  the  numerical 
experiments  that  the  use  of  periodic  extension  gives  accurate  results.  In  addition,  it 
allows  decoupling  of  the  ODEs  using  eigenvalue  analysis  and  thus  reduces  the 
computational  cost.  Here,  after  expressing  the  unknown  coefficients  lying  outside  the 
finite  domain  in  terms  of  the  inner  coefficients  considering  periodic  extension,  the  ODEs 
given  by  Eq.  (21)  can  be  written  as  a  matrix  equation  of  the  form 
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where  A*  is  the  first-order  connection  coefficient  matrix  obtained  after  periodic 
extension.  The  coupled  ODEs  given  by  Eq.(22)  are  decoupled  using  eigenvalue  analysis 
similar  to  that  done  in  time  approximation.  It  should  be  mentioned  here  that  matrix  A' 
obtained  after  periodic  extension  has  a  circulant  form  and  its  eigen  parameters  are  known 
analytically  [10].  Let  the  eigenvalues  be  jy5,  ,then  the  decoupled  ODEs  corresponding  to 
Eq.  (22)  are 
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where  and  similarly  other  transformed  displacements  are 

iToj  =  (24) 

where  T*  is  the  eigenvector  matrix  of  A* . 


Following  exactly  similar  steps,  the  final  transformed  and  decoupled  form  of  the  Eqs.  (7) 
and  (8)  (following  reduction  using  temporal  approximation)  are 
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Similarly,  the  transformed  form  of  the  force  boundary  conditions  given  by  Eqs.  (9)-(12) 
(following  reduction  using  temporal  approximation)  are 
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The  final  transformed  ODEs  given  by  Eqs.  (23), (25),  and  (26)  and  the  boundary 
conditions  Eqs.  (28)-(30)  are  used  for  2D  WSEM  formulation  and  is  similar  to  the  2D 
ESEE  technique  [1]. 

Similar  to  the  temporal  approximation,  the  spatial  sampling  rate  AT  is  also  determined 
from  the  order  N  of  the  scaling  function  used  for  spatial  approximation  and  the  spatial 
distribution  of  the  load. 


The  four  degrees  of  freedom  per  node  associated  with  the  element  formulation  are  , 
Vgij ,  iv^  and  /  dx  as  shown  in  Eigl. 


Figure  1:  Spectral  Plate  Element  with  all  the  degrees  of  freedom  and  Stress 

Resultants 


The  corresponding  nodal  forces  are  ,  M  and  Vy .  From  the  previous  sections, 

for  unrestrained  lateral  edges,  we  get  a  set  of  decoupled  ODEs  (Eqs.  (23),  (25),  and  (26)) 
for  an  isotopic  plate  using  CLPT,  in  a  transformed  wavelet  domain.  These  equations  are 
required  to  be  solved  for  u^y  ,  v^y ,  w^y ,  and  the  actual  solutions  uo(x,  y  ,  t),  vo(X,  y  ,  t), 

w(x,  y  ,  t)  are  obtained  using  inverse  wavelet  transform  twice  for  spatial  Y  dimension  and 
time. 

It  can  be  seen  that  the  transformed  decoupled  ODEs  have  a  form  that  is  similar  to  that  in 
FSEM[1],  and  thus,  the  formulation  of  WSEM  from  here  is  similar  to  FSEM  formulation 
or  1-D  WSEM  formulation  given  in  Reference  [5]  .  Thus,  the  formulation  is  not  repeated 

here.  Finally,  the  transformed  nodal  forces  and  transformed  nodal  displacements 
{m  }  are  related  as 


{r}  =  [K^]{u"}  (31) 

where  is  the  exact  elemental  dynamic  stiffness  matrix.  The  solution  of  the  Eq.  (31) 

and  the  assembly  of  the  elemental  stiffness  matrices  to  obtain  the  global  stiffness  matrix 
are  exactly  similar  to  conventional  FE  technique. 

4.  NUMERICAL  EXAMPLES. 

Here,  the  axial  and  transverse  wave  propagations  in  composite  graphite-epoxy  plates  of 
different  configurations  and  ply  orientations  are  studied.  The  examples  used  (shown  in 
Figs.  2  (a),  (b)  and  (c))  consist  of  uniform  cantilever  plate  with  symmetric  and 
asymmetric  ply  lay-ups,  cantilever  asymmetric  plate  with  ply  drop  and  symmetric  folded 
plate  structure.  The  material  properties  are  as  follows,  El  =  144.48  GPa,  E2  =  E3  =  9.63 
GPa,  G23  =  G13  =  G12  =  4.128  GPa,  v23  =  0.3,  vl3  =  vl2  =  0.02,  p  =  1389  kg/ml 


Figure  2:  (a)  Uniform  cantilever,  (b)  Ply  drop  and  (c)  Folded  plates. 

In  all  the  examples  provided,  the  load  applied  is  an  unit  impulse  and  is  shown  in  time  and 
frequency  domains  in  Fig.3.  The  load  is  applied  at  the  edge  along  the  Y  -axis  and  has  a 
pulse-like  spatial  distribution  and  is  given  as 

FiY)  = 

where,  a  is  a  constant  and  can  be  varied  to  change  the  Y  -axis  variation  of  the  load. 


Time  (^ls) 


Figure  3:  Impact  load  and  Fourier  transform  of  the  load  (inset) 


The  2-D  WSEM  model  is  formulated  with  Daubechies  scaling  function  of  order  N  =22 
for  temporal  approximation  and  =  4  for  spatial  approximation.  The  time  sampling  rate 
is  At  =  2  jus,  unless  otherwise  mentioned,  while  the  spatial  sampling  rate  zIT  is  varied 
depending  on  LY  and  load  distribution  F(Y ). 

4.1  Uniform  Plate 

The  uniform  cantilever  plate  is  shown  in  Fig.  2(a)  and  is  fixed  at  one  edge  CD  and  free  at 
the  other  edge  AB  along  Y  -axis.  Numerical  experiments  are  performed  by  considering 
the  other  two  edges  AC  and  BD  along  X-axis  to  be  free-free.  The  dimensions  are  LX  and 
LY  along  X  and  Y  axes  respectively,  while  the  depth  (=  2h)  is  kept  fixed  at  0.01  m  with  8 
laminates.  The  spectrum  relations  for  the  plate  with  LY  =  0.25  m  and  symmetric  ply 
orientation  of  [02/904/02]  are  plotted  in  Fig.  4  .  Figs.  4  (a)  and  (b)  respectively  show  the 
real  and  imaginary  parts  of  the  wavenumbers  for  a  Y  -wavenumber  of  50.  It  can  be  seen 
that  the  wavenumber  has  significant  real  and  imaginary  parts.  This  implies  that  the  waves 
are  inhomogeneous  in  nature,  i.e  they  attenuate  as  they  propagate.  The  wavenumbers 
have  been  obtained  with  Jt  =  4  //s  i.e  for  a  Nyquist  frequency  of  Jnyq  =  125  kHz.  As  said 
earlier,  WSEM  predicts  accurate  wavenumbers  only  up  to  a  certain  fraction  pN  of 
Nyquist  frequency  This  fraction  pN  for  A  =  22  is  approximately  0.6.  Thus  in  Figs. 
4. (a)  and  (b),  the  wavenumbers  are  plotted  up  to  a  frequency  fN  =  pNfnyq  =  75  kHz.  In 
Figs.  5(a)  and  (b),  the  real  and  imaginary  parts  of  the  wavenumber  of  the  composite  plate 
with  LY  =  0.25  m  and  asymmetric  ply  lay  up  of  [04/904]  are  presented  respectively.  The 
plots  are  obtained  for  a  Y  wavenumber  50  and  time  sampling  rate  At  =  A  ps  or  Nyquist 
frequency  of fnyq  =  125  kHz.  Even  here,  the  wavenumbers  calculated  using  the  present 
method  are  correct  up  to  a  frequency  of  fN  =  pNfnyq  =  75  kHz.  Similar  to  plate  with 
symmetric  ply  lay  up,  the  wavenumbers  for  asymmetric  ply  orientation  have  considerable 
real  and  imaginary  parts  and  thus  the  waves  are  inhomogeneous  in  nature.  It  can  be  seen 
that  for  both  symmetric  and  asymmetric  plates,  the  real  part  of  the  wavenumbers  show 
three  propagating  modes  corresponding  to  axial,  lateral  and  transverse  wave  propagation. 
There  are  three  cut-off  frequencies  which  vary  with  wavenumbers. 


Figure  4:  The  (a)  real  and  (b)  imaginary  parts  of  the  wavenumber  of  plate  with 

symmetric  ply  lay-up  of  [02/904/02]. 


Figure  5:  The  (a)  real  and  (b)  imaginary  parts  of  the  wavenumber  of  plate  with 

asymmetric  ply  lay-up  of  [04/904]. 

Next,  the  time  domain  responses  of  a  plate  with  LX  =  0.5  m,  LY  =  0.25m  and  a  symmetric 
ply  lay  up  of  [08],  simulated  using  WSFE  method  are  validated  with  2-D  FE  results.  In 
Figs.  6  (a)  and  (b),  respectively,  the  axial  and  transverse  velocities  measured  at  the  mid 
point  of  edge  AB  (see  Fig.  2  (a))  of  the  cantilever  plate  are  plotted  and  compared  with  2- 
D  FE  results.  The  impulse  load,  shown  in  Fig.  3,  is  applied  along  AB,  correspondingly  in 
axial  and  transverse  directions.  The  Y  variation  of  the  load  is  obtained  using  a  =  0.03.  As 
mentioned  earlier  only  one  WSEM  is  used  to  model  the  structure  and  the  time  window  is 
kept  Tw  =  5\2  //s  with  number  of  sampling  points  n  =  256  and  At  =  2  jus.  The  number  of 
discretization  points  along  Y  -axis  is  m  =  64  and  thus  the  spatial  sampling  rate  is  AY  =  LY 
/(m  -  1)  =  0.004  m. 


Figure  6:  (a)  Axial  and  (b)  transverse  velocities  at  mid-point  of  edge  AB  in  a  [0]8 
cantilever  plate  (see  Fig.  2(a))  with  LX  =  0.5  m  and  LY  =  0.25  m  due  to  tip  impulse 
load  applied  in  axial  and  transverse  directions  along  AB  respectively 

A  very  refined  mesh  with  6432,  4-noded  plane  stress  quadrilateral  elements  is  used  for 
the  2-D  FE  analysis,  while,  Newmark’s  scheme  with  time  step  1  ^s  was  used  for  time 


integration.  It  can  be  seen  that  WSEM  and  FE  results  match  very  well.  A  comparison  is 
also  provided  with  FSEM  results.  As  stated  earlier,  it  can  be  seen  from  the  figures,  that 
unlike  WSEM,  FSEM  is  unable  to  accurately  capture  the  reflections  from  the  lateral 
edges  AC  and  BD  in  this  example.  The  velocities  obtained  from  FSEM  modeling  show 
only  the  reflections  from  the  fixed  edge  CD.Thus,  for  structures  with  finite  or  short 
dimensions,  FSEM  results  will  deviate  substantially  from  the  actual  responses.  In 
addition,  simulation  with  FSFE  requires  ''throw-ojf'  element  to  impart  artificial  damping 
to  the  structure  and  a  large  time  window  Tw  =  16384  ^s  {At  =  2  jus  and  n  =  8192)  to 
remove  the  distortions  due  to  ''wrap  around"  problem.  It  should  be  restated  here,  that  the 
accuracy  of  the  response  simulated  using  WSEM  is  independent  of  the  time  window  Tw 
which  is  chosen  as  required  for  observation.  Similar  validation  with  2-D  FE  results  are 
done  for  the  plate  with  asymmetric  ply  orientation  of  [04/904]  with  the  other 
configurations  and  loading  condition  remaining  same  as  the  previous  example.  The 
parameters  for  WSEM  modeling  and  the  FE  mesh  are  similar  to  that  used  in  the  last 
example  for  symmetric  lay-up  plate.  Figs.  7(a)  and  (b)  show  the  axial  and  transverse 
velocities  at  mid-point  of  the  free  edge  AB  due  loading  in  axial  and  transverse  directions 
respectively.  In  this  case  also  the  responses  obtained  with  WSFE  match  well  with  FE 
results.  In  addition,  the  corresponding  velocities  simulated  using  FSEM  using  the  same 
time  window  Tw,  time  At  and  spatial  AY  rates,  are  also  plotted  to  show  the  limitations  of 
the  FSEM  method  in  modeling  finite  dimension  2-D  structures. 


Figure  7:  (a)  Axial  and  (b)  transverse  velocities  at  mid-point  of  edge  AB  in  a 
[04/904]  cantilever  plate  (see  Fig.  4.11(a))  with  LX  =  0.5  m  and  LY  =  0.25  m  due  to 
tip  impulse  load  applied  in  axial  and  transverse  directions  along  AB  respectively. 


Figure  8:  (a)  Axial  and  (b)  transverse  velocities  at  mid-point  of  edge  AB  in 
[02/904/02],  [02/604/02]  and  [0]8  cantilever  plates  (see  Fig.  4.11(a))  with  LX  =  0.5  m 
and  LY  =  0.25  m  due  to  tip  impulse  load  applied  in  axial  and  transverse  directions 

along  AB  respectively. 


In  Fig  8(a),  the  axial  velocities  at  the  mid  point  of  edge  AB  of  cantilever  plate  shown  in 
Fig.  2(a)  are  plotted  for  two  different  symmetric  ply  orientations  of  [02/904/02]  and 
[02/604/02],  The  plate  has  a  dimension  of  LX  =  0.5  m  and  LY  =  0.25  m.  The  impulse 
loading  is  similar  to  that  in  previous  example  and  is  applied  in  axial  direction  along  edge 
AB  with  a  =  0.03  for  obtaining  the  Y  variation  of  load.  Even  here,  the  simulations  are 
done  with  a  single  WSEM,  time  sampling  rate  At  =  2  //s,  time  window  Tw  =  512  //s,  m  = 
64  and  thus  Y  sampling  rate  dT  =  0.004  m.  The  corresponding  transverse  velocities  at  the 
mid  point  of  edge  AB  due  transverse  loading  are  plotted  in  Fig.  8(b).  It  can  be  seen  that 
the  amplitude  of  the  incident  axial  waves  is  more  for  the  ply  lay  up  [02/904/02]  than  that 
for  [02/604/02].  This  is  expected  as  the  longitudinal  stiffness  of  the  lay  up  [02/904/02]  is 
more  than  [02/904/02].  However,  the  amplitudes  do  not  vary  much  for  the  transverse 
wave  velocities.  For  similar  reason  the  reflections  from  the  lateral  edges  appear  later  in 
[02/904/02]  than  in  [02/604/02].  The  axial  velocities  at  the  mid  point  of  edge  AB  of  the 
cantilever  plate  shown  in  Fig.  2(a)  with  asymmetric  ply  lay-ups  [04/904]  and  [04/604]  are 
plotted  in  Fig.  9(a). 


Figure  9:  (a)  Axial  and  (b)  transverse  velocities  at  mid-point  of  edge  AB  in  [04/904] 
and  [04/604]  cantilever  plates  (see  Fig.  4.11(a))  with  LX  =  0.5  m  and  LY  =  0.25  m 
due  to  tip  impulse  load  applied  in  transverse  and  axial  directions  along  AB 

respectively. 

These  responses  are  due  to  impulse  load  (see  Fig.  3)  applied  along  AB  in  transverse 
direction  and  result  from  the  axial-flexural  coupling  arising  due  to  the  asymmetric  ply 
la3mp.  The  Y  variation  of  the  load  is  obtained  with  cc  =  0.03.  The  plate  dimensions  are  LX 
=  0.5  m  and  LY  =  0.25  m.  The  parameters  involved  in  WSEM  modeling  are  same  as  in 
previous  example.  Similar  to  the  symmetric  configurations,  the  amplitude  of  incident 
wave  for  [04/904]  is  more  than  [04/604]  due  to  reasons  as  explained  before.  In  Fig.  9(b) 
the  transverse  velocities  at  mid-point  of  edge  AB  due  to  the  load  applied  in  axial 
direction  are  plotted.  These  responses  also  result  from  the  axial-flexural  coupling. 
However,  the  amplitudes  of  these  velocities  arising  due  to  the  coupling  vary  considerably 
with  the  ply  orientations  and  as  expected  are  higher  for  [04/904]  due  to  its  lower  stiffness. 

Figs.  10  (a)  and  (b)  show  the  snapshots  of  the  axial  velocities  of  the  cantilever  plate 
shown  in  Fig.  2(a)  with  a  symmetric  ply  orientation  of  [08],  at  time  instances  T  =  250  //s 
and  T  =  375  fis  respectively.  The  plate  dimensions  are  LX  =  2.0  m  and  LY  =  0.5  m,  and  is 
modeled  using  single  WSEM  with  m  =  64  sampling  points  in  Y  direction.  The  impulse 
load  as  in  Eig.  3  is  applied  along  edge  AB  in  axial  direction  and  the  Y  variation  is 
obtained  with  a  =  0.05.  The  snapshot  at  T  =  250  //s  shows  the  forward  moving  axial 
wave,  while  at  T  =  375  //s  the  wave  reflected  from  the  fixed  end  of  the  cantilever  plate 
can  be  captured.  It  should  be  mentioned  here  that  the  velocities  at  all  the  sampling  points 
along  Y  direction  and  at  any  point  along  X  direction  used  to  obtain  the  snapshots  are 
obtained  from  a  single  simulation.  In  Figs.  11(a)  and  (b)  the  snapshots  of  the  transverse 
velocities  of  the  cantilever  plate  shown  in  Fig.  2(a)  with  a  symmetric  ply  orientation  of 
[08]  are  presented  for  time  instances  T  =  500  /is  and  T  =  1000  //s  respectively.  The 
dimensions  and  the  loading  conditions  are  kept  same  except  that  the  loading  here  is 
applied  in  transverse  direction. 


(a)  (b) 


Figure  10  :  Snapshots  of  axial  velocities  at  time  instances  (a)  T  =  250  //s  and  (b)  T  = 
375  fis  in  a  [0]8  cantilever  plate  (see  Fig.  2(a))  with  LX  =  2.0  m  and  LY  =  0.5  m  due 
to  tip  impulse  load  applied  in  axial  direction. 


(a)  (b) 

Figure  11:  Snapshots  of  transverse  velocities  at  time  instances  (a)  T  =  500  //s  and  (b) 
T  =  1000  //s  in  a  [0]8  cantilever  plate  (see  Fig.  2(a))  with  LX  =  2.0  m  and  LY  =  0.5  m 
due  to  tip  impulse  load  applied  in  transverse  direction 

WSEM  modeling  is  also  performed  with  the  similar  parameters.  Here,  both  the  snapshots 
show  the  forward  propagation  of  the  flexural  waves  and  the  dispersive  nature  of  these 
waves  can  also  be  observed. 

4.2  Ply  Dropped  and  Folded  Plates  Analysis 

Ply  dropping  is  done  to  taper  a  composite  plate  and  is  in  very  much  use  in  many 
composite  structures.  This  example  of  wave  propagation  in  ply  dropped  plate  due  to 
impulse  load  is  performed  to  emphasize  the  effectiveness  of  the  present  method  in 
modeling  such  complex  structures.  The  configuration  of  the  plate  is  shown  in  Fig.  2(b) 
and  the  overall  dimension  is  LX  =  1.5  m  in  X  direction  and  LY  =  0.5  m  in  F  direction.  It 
is  fixed  at  one  end  CD  and  free  at  the  other  edge  AB.  The  two  lateral  edges  AC  and  BD, 
parallel  to  Y  axis  are  considered  free.  The  plate  is  divided  into  three  regions  along  X 
direction.  From  the  fixed  end  to  0.5  m,  the  ply  lay-up  is  [04/904]  and  the  total  thickness  is 
h\  =0.01  m.  For  the  next  0.5  m,  two  plies  are  reduced.  Here  the  thickness  is  hi  =  0.0075 
m  and  the  lay-up  is  [03/903].  The  last  part  has  a  thickness  of  h3  =  0.005  m  and  the  lay-up 


is  [02/902],  As  there  are  two  discontinuities  present,  three  WSEM  are  assembled  to 
model  the  plate  and  the  number  of  sampling  points  along  Y  direction  is  m  =  64  with  JF  = 
0.008  m.  The  unit  impulse  load  (see  Fig.  3)  is  applied  along  edge  AB  and  the  Y  variation 
of  the  load  is  obtained  with  a  =  0.05.  In  Fig.l2  (a),  the  axial  velocity  measured  at  the 
mid-point  of  edge  AB  due  to  load  in  axial  direction  is  plotted. 


Figure  12:  (a)  Axial  and  (b)  transverse  velocities  of  asymmetric  [04/904]  uniform 
and  ply  dropped  plate  (Fig.  2(b)  and  (c)  respectively)  with  LX  =  1.5  m,  LY  =  0.5  m 

along  edgeAB. 

The  axial  velocity  of  corresponding  uniform  cantilever  plate  without  ply  drop  and  LX  = 
1.5  m  and  LY  =  0.5  m,  is  also  plotted  for  comparison.  As  expected,  the  amplitude  of  the 
incident  wave  for  the  uniform  plate  is  considerably  less  than  that  for  ply  dropped  plate.  In 
addition,  there  are  also  much  differences  in  the  reflected  waves  arriving  later.  This  is 
because,  apart  from  the  reflections  from  the  lateral  free  edges  and  the  fixed  edge,  the 
response  of  the  ply  dropped  plate  consists  of  reflections  from  the  discontinuities  present. 
These  reflections  are  absent  in  the  response  of  the  uniform  plate  In  Fig.  12(b),  the 
corresponding  transverse  velocities  due  to  transverse  impulse  load,  are  plotted  for  both 
ply  dropped  and  uniform  plates.  Similar  to  the  axial  velocities,  here  also  there  is  a 
prominent  difference  in  the  amplitudes  of  the  incident  waves  for  the  uniform  and  ply 
dropped  plates.  As  explained  earlier,  the  reflected  waves  also  show  considerable 
variations. 

Next,  wave  propagation  analysis  in  a  folded  plate  structure,  shown  in  Fig.  2(c),  is  done 
using  the  formulated  WSFM  method.  The  folded  plate  model  is  obtained  by  assembling 
three  2-D  WSFM.  The  number  of  sampling  point  along  y  direction  is  m  =  32  for  all  the 
three  plates.  The  plates  have  the  same  dimensions  of  Lx  =0.5  m  and  Ly  =  0.5  m  and  a 
symmetric  ply  lay-up  of  [08].  The  two  slanted  plates  are  at  an  angle  of  60“  to  the  X 
direction  and  the  two  ends  denoted  as  AG  and  DH  are  considered  fixed.  Assembling  is 
done  along  the  local  x  direction  (see  Fig.  2(c))  after  the  required  transformation  similar  to 
conventional  1-D  FF.  The  unit  impulse  load  shown  in  Fig.  3  is  applied  along  edge  BF 
and  CF.  The  Y  variation  of  the  load  is  obtained  with  a  =  0. 1 . 


Figure  13:  (a)  Axial  and  (b)  transverse  velocities  at  mid  point  of  BC  of  [08]  folded 
plate  (Fig.  2(c))  due  load  applied  along  edge  BE  and  CF  in  axial  and  transverse 

directions  respectively. 

In  Figs.  13  (a)  and  (b),  the  axial  and  transverse  velocities  measured  at  the  mid-point  of 
the  edge  BC,  due  to  the  above  mentioned  load  applied  in  axial  (X)  and  transverse  (Z) 
directions  are  plotted  respectively. 

These  numerical  experiments  show  the  capability  of  the  present  method  to  analyze  the 
rather  complicated  folded  plate  structures.  It  should  be  mentioned  here  that  2-D  FSEM 
cannot  be  used  for  modeling  such  finite  dimension  structures.  In  addition  even  for  semi 
infinite  structures,  the  responses  simulated  with  2-D  FSEM  will  be  highly  distorted  due  to 
''wrap  around'  as  there  are  multiple  reflections  from  the  several  discontinuities  present 

5.  CONCLUSIONS 

The  main  aim  of  this  work  is  to  bring  out  the  wave  propagation  behavior  in  composite 
plates.  First  the  wavelet  based  spectral  finite  element  is  developed  for  composite  plate 
with  different  ply  orientations.  The  developed  WSEM  is  used  to  study  both  frequency 
and  time  domain  wave  characteristics  such  as  the  spectrum  and  dispersion  relations.  The 
time  domain  responses  are  compared  with  2-D  FE  analysis  results  and  they  match  well. 
Comparisons  are  also  provided  with  the  2-D  FSEM  analysis  results  to  highlight  the 
efficiency  of  2-D  WSEM  for  analysis  of  finite  dimension  plates  unlike  FSEM.  Apart 
from  uniform  cantilever  plate,  responses  are  also  simulated  for  plates  with  ply-drop  and 
folded  structures.  These  examples  show  the  effectiveness  of  the  proposed  methodology  in 
modeling  rather  complex  structures.  The  next  logical  extension  of  this  work  is  to  model 
the  damages  such  as  delaminations  and  fibre  breakages  in  composites.  The  delamination 
can  be  an  area  delamination  or  through  the  depth  delamination.  First  the  latter  will  be 
attempted  in  the  next  year  of  the  project.  Similarly,  through  the  width  transverse  cracks 
simulating  the  fiber  breakage  will  also  be  undertaken  in  coming  year. 
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